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Abstract 

At  HF  through  UHF  frequencies,  wave  propagation  in  a  forest  environment 
is  mainly  attributed  to  a  lateral  wave  which  propagates  at  the  canopy-air  in¬ 
terface.  Due  to  the  existence  of  tree  trunks,  significant  multiple  scattering  also 
occurs  which  is  the  dominant  source  of  field  fluctuations.  Basically  the  cur¬ 
rent  induced  in  the  tree  trunks  by  the  source  and  the  lateral  wave  re-radiate 
and  generate  higher  order  lateral  waves  and  direct  scattered  waves.  Using  a 
full-wave  analysis  based  on  the  method  of  moments  in  conjunction  with  Monte 
Carlo  simulations,  the  effect  of  multiple  scattering  among  a  very  large  number 
of  tree  trunks  is  studied.  It  is  shown  that  only  scatterers  near  the  source  and 
the  observation  points  contribute  to  the  field  fluctuations  significantly.  This 
result  drastically  simplifies  the  numerical  complexity  of  the  problem.  Keeping 
about  200  tree  trunks  in  the  vicinity  of  the  transmitter  dipole  and  the  receiver 
point,  a  Monte  Carlo  simulation  is  used  to  evaluate  the  statistics  of  the  spatial 
and  spectral  behavior  of  the  field  at  the  receiver.  Using  a  wideband  simulation, 
the  temporal  behavior  (impulse  response)  is  also  studied  as  is  performance  of 
antenna  arrays  and  the  effects  of  different  spatial  diversity  combining  schemes 
in  such  a  multi-path  environment. 


1  Introduction 

Accurate  prediction  of  radio  wave  propagation  in  a  communications  channel  is  es¬ 
sential  in  the  development  and  design  of  an  efficient,  and  low-cost  wireless  system. 
High  fidelity  models  are  necessary  in  order  to  treat  the  tradeoff  between  radiated 
power  and  signal  processing  by  addressing  such  issues  as  coherency,  field  variations, 
multi-path,  and  dispersive  (path  delay)  effects.  Current  channel  models  are  heuristic 
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in  nature  and  have  limited  applicability.  A  physics-based  approach  to  channel  charac¬ 
terization  gives  insight  into  the  mechanisms  of  radio  wave  propagation  and  inherently 
provides  highly  accurate  results.  With  this  as  a  motivation,  a  physics-based  model¬ 
ing  approach  is  being  pursued  at  the  University  of  Michigan  for  which  advanced  and 
efficient  electromagnetic  diffraction  models  are  being  developed. 

Due  to  relatively  high  attenuation  rates,  direct  wave  propagation  in  a  forest  envi¬ 
ronment  is  not  possible  over  large  distances  at  high  frequencies.  In  the  HF-UHF  range 
where  both  the  transmitter  and  receiver  are  embedded  in  the  foliage,  radio  signals 
can  propagate  over  relatively  large  distances.  This  peculiar  behavior  is  explained  by 
certain  type  of  surface  waves  known  as  the  lateral  wave  [1]. 

Mathematically  speaking,  the  lateral  wave  is  the  contribution  of  the  branch  cut 
from  the  integral  of  the  spectral  representation  of  a  dipole  field  inside  a  half-space 
dielectric.  This  formulation  provides  an  expression  for  the  mean-field  assuming  the 
canopy-air  interface  is  smooth.  For  predicting  the  mean-field  more  accurately,  in  a 
recent  study  [2],  the  effect  of  canopy-air  interface  roughness  on  the  propagation  of 
lateral  waves  was  studied  and  it  was  shown  that  the  canopy-air  interface  roughness 
reduces  the  mean-field.  As  mentioned  earlier  path-loss  only  partially  characterizes 
the  channel  and  other  scattering  mechanisms  must  be  accounted  for  to  complete  the 
model.  At  UHF  and  lower  frequencies,  the  dimensions  of  leaves  and  branches  are 
small  compared  to  wavelength  and  do  not  cause  significant  scattering  which  is  the 
source  of  signal  fluctuations,  multi-path,  and  dispersion.  The  effect  of  tree  trunks 
which  are  electrically  large  create  a  highly  scattering  environment. 

In  this  paper  the  effect  of  tree  trunks  is  accounted  for  by  computing  their  interac¬ 
tion  with  the  source  and  the  lateral  waves.  A  numerical  solution  based  on  the  method 
of  moments  (MoM)  in  conjunction  with  Monte-Carlo  simulation  is  proposed  to  eval¬ 
uate  the  scattering  effects  of  tree  trunks.  However,  considering  the  number  of  tree 
trunks  between  the  transmitter  and  receiver,  it  is  quite  obvious  that  a  brute-force  ap¬ 
plication  of  MoM  is  not  possible  due  to  the  exorbitant  memory  and  computation  time 
requirements.  To  make  the  solution  tractable  while  maintaining  the  model  fidelity, 
three  techniques  are  proposed:  1)  simplification  of  the  MoM  formulation  noting  that 
tree  trunks  are  sparsely  distributed,  2)  simplification  base  on  physical  insight  by  not¬ 
ing  that  the  scattered  fields  from  tree  trunks  between  the  transmitter  and  receiver, 
but  distant  from  them,  are  almost  in-phase,  and  3)  simplification  of  field  computation 
using  the  reciprocity  theorem.  In  what  follows,  the  forest  model  is  described  first  and 
then  the  simplified  MoM  formulation  is  presented.  This  model  is  used  to  demonstrate 
that  only  scatterers  near  the  source  and  observation  points  significantly  contribute  to 
the  field  fluctuations.  Next  the  near-field  interaction  of  tree  trunks  with  the  source 
field  and  lateral  waves  are  computed  by  application  of  reciprocity  and  using  the  MoM 
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formulation.  Finally,  numerical  simulations  are  presented,  where  the  spatial  decor¬ 
relations  in  a  forest  environment  are  examined  and  the  performance  of  the  antenna 
arrays  and  spatial  diversity  schemes  evaluated. 

2  Forest  Model 

The  model  presented  in  this  section  is  suitable  for  predicting  the  statistics  of  the  field 
radiated  by  an  elementary  antenna  embedded  in  a  forest  environment  and  is  valid 
for  frequencies  up  to  UHF.  As  mentioned  earlier,  since  the  dimensions  of  leaves  and 
branches  are  small  compared  to  the  wavelength,  the  forest  canopy  can  be  modeled  by 
a  homogeneous  dielectric  .  However,  the  trunks  whose  dimensions  are  comparable  to 
or  larger  than  a  wavelength  must  be  treated  separately.  Figure  1  shows  the  geometry 
of  the  wave  propagation  problem  where  both  the  source  and  observation  points  are 
within  a  dielectric  slab  with  effective  permittivity  £ejf-  This  dielectric  slab  is  assumed 
to  have  a  smooth  lower  interface  with  a  dielectric  half-space  (representing  the  ground) 
and  a  rough  upper  interface  with  air.  Tree  trunks  are  modeled  by  dielectric  cylinders 
perpendicular  to  the  ground  plane  having  a  relatively  large  height-to-diameter  ratio. 
Ignoring  the  scattering  from  tree  trunks,  the  mean  field  at  the  receiver  is  composed 
of  the  following  components  as  depicted,  in  Fig.  2  : 

1.  geometric  optics  terms  which  include  the  direct  propagation  between  the  trans¬ 
mitter  and  receiver  and  reflected  fields  from  the  upper  and  lower  interfaces.  The 
mean  reflected  field  from  the  upper  interface  is  reduced  exponentially  where  the 
exponent  is  proportional  to  the  rms  height  of  the  canopy-air  interface  roughness 
[3].  These  terms  are  only  important  when  the  receiver  is  relatively  close  to  the 
transmitter,  otherwise  due  to  the  lossy  nature  of  the  effective  dielectric  constant 
of  the  foliage  these  components  do  not  contribute  much.  Basically  the  geomet¬ 
ric  optics  terms  exponentially  decay  with  distance  between  the  transmitter  and 
the  receiver. 

2.  Lateral  wave  (a  diffraction  term)  which  propagates  along  the  canopy-air  inter¬ 
face  and  decays  proportionally  to  the  reciprocal  of  the  radial  distance  squared (~ 
1/p2).  The  path  loss  associated  with  the  lateral  wave  propagation  increases  with 
increasing  foliage  density  (effective  dielectric  constant)  and  decreases  by  bring¬ 
ing  the  source  or  observation  points  closer  to  the  canopy-air  interface.  The 
path  loss  increases  as  the  canopy-air  interface  roughness  (rms  height)  relative 
to  wavelength  increases.  Close  examination  of  the  expression  representing  the 
lateral  wave  contribution  reveals  that  the  contribution  can  be  attributed  to  a 
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ray  that  emanates  from  the  source  along  the  critical  angle,  propagates  along  the 
canopy-air  interface,  and  arrives  at  the  receiver  along  the  critical  angle  as  shown 
in  Fig.  2.  Other  higher  order  lateral  wave  contributions  also  exist,  of  which 
two  are  significant.  One  corresponds  to  the  ray  emitted  from  the  source  which 
reflects  from  the  ground  boundary,  propagates  along  the  critical  angle  and  then 
travels  along  the  interface  ( Eql )•  The  other  one  is  the  reflected  lateral  wave 
from  the  ground  which  arrives  at  the  receiver  ( Elg )•  These  two  contributions 
are  also  shown  in  Fig.  2 


The  expression  for  the  electric  field  resulted  from  the  asymptotic  expansion  of  the 
spectral  integral  around  the  branch  cut  (Lateral  wave)  for  a  half-space  dielectric  with 
a  smooth  boundary  is  given  by  [2] 
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where  p  is  the  radial  distance  between  the  transmitter  and  receiver,  h  and  h!  are  the 
depth  of  transmitter  and  receiver  below  the  canopy-air  interface,  II  is  the  current 
moment  of  the  dipole,  and  k0  and  Zq  are,  respectively,  the  propagation  constant  and 
the  characteristic  impedance  of  free-space.  In  (1)  the  unit  vector,  /,  denotes  the  dipole 
orientation  k  —  l/eeff  where  £efj  is  the  effective  dielectric  constant  of  foliage,  and  A 
is  a  symmetric  dyad  given  by 
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Here  <j>  is  a  cylindrical  coordinate  angle  representing  the  location  of  observation  point 
with  respect  to  the  transmitter.  Assuming  the  canopy  height  is  H,  the  lateral  wave 
from  the  image  of  the  source  in  the  ground  plane  ( Eql )  is  given  by 
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and  R±  and  are  the  Fresnel  reflection  coefficients  of  the  ground  evaluated  at  the 
critical  angle  6C  =  sin-1  s/R.  Similarly  ELG  is  given  by 
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In  (3)  and  (5)  the  effects  of  the  ground  surface  wave  is  ignored.  The  phase  term  of 
equation  (1)  indicates  that  the  lateral  wave  is  locally  a  plane  wave  propagating  along 


the  direction  kx  =  k0 


and  also  noting  EL  ■  kx  =  0.  The 


|cos  (j)x  +  sin  (fry  — 

expression  for  the  lateral  wave  contribution  for  rough  canopy-air  interface  is  more 
complicated  than  (1),  but  the  mean-field  still  represents  a  plane  wave  locally  [2].  To 
include  the  effects  of  scattering  from  tree  trunks  consider  the  geometry  of  the  problem 
as  shown  in  Fig.  3.  The  source  and  its  image  excite  polarization  currents  in  the  di¬ 
electric  cylinders  which  in  turn  re-radiate  and  produce  multiple  scattering  among  the 
tree  trunks  and  secondary  lateral  waves  (lateral  waves  generated  by  scattering  from 
the  tree  trunks)  that  arrive  at  the  receiver.  The  sum  of  all  lateral  and  the  secondary 
lateral  waves  will  be  referred  to  as  the  total  incident  wave  in  the  vicinity  of  all  receiver. 
The  total  incident  wave  excites  polarization  currents  within  the  dielectric  cylinders 
(tree  trunks)  near  the  receiver,  which  re-radiate  and  together  with  the  total  incident 
fields  constitute  the  total  field  at  the  receiver.  A  formal  solution  to  the  problem  can 
be  obtained  rather  easily  by  casting  the  formulation  in  terms  of  an  integral  equa¬ 
tion  for  the  polarization  currents  induced  inside  the  dielectric  cylinders.  However, 
the  brute-force  solution  of  the  integral  equation  using  the  method  of  moments  is  not 
tractable  because  of  the  large  number  of  unknowns  and  the  complex  nature  of  the 
dyadic  Green’s  function  of  the  problem.  To  make  the  problem  tractable  an  accurate 
approximate  solution  is  sought.  An  accurate  approximate  solution,  considering  the 
physics  of  the  problem,  can  be  obtained  as  will  be  shown  in  the  following  section. 


3  Reduced  Problem 

As  mentioned  before,  estimation  of  field  fluctuations  in  a  forest  environment  requires 
the  computation  of  multiple  scattering  effects,  as  well  as  the  interaction  of  lateral 
waves  with  large  number  of  dielectric  cylinders.  It  is  expected  that  the  induced 
polarization  currents  in  cylinders  near  the  source  be  much  stronger  than  those  in 
cylinders  distant  from  the  source.  Also  the  contribution  to  the  received  fields  from 
cylinders  in  the  vicinity  of  the  receiver  is  expected  to  be  high.  Although  the  contri¬ 
bution  to  the  scattered  fields  from  individual  cylinders  between  the  transmitter  and 
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receiver  which  are  not  in  the  close  proximity  of  the  either  is  relatively  small,  one 
may  argue  that  there  are  many  such  cylinders  and  the  overall  contribution  may  be 
significant.  Considering  the  path-length  between  the  transmitter  to  a  cylinder  and 
from  the  cylinder  to  the  receiver,  and  noting  that  the  scattering  is  strongest  in  near 
forward  direction,  the  scattered  fields  from  these  cylinders  arrive  almost  in-phase. 
Therefore  the  scattered  field  from  cylinders  that  are  not  close  to  the  receiver  and 
transmitter  are  not  expected  to  contribute  to  the  field  fluctuations  significantly.  The 
effect  of  these  scatterers  may  be  accounted  for  by  replacing  them  with  an  effective 
dielectric  constant.  To  examine  the  validity  of  these  postulations,  a  2-D  medium  con¬ 
sisting  of  infinite  cylinders,  excited  by  line  source,  is  considered.  Figure  4(a)  shows 
the  geometry  of  the  problem  where  a  line  source  capable  of  supporting  z-directed 
current  (TM  case)  or  p-directed  current  (TE  case)  is  used  as  the  transmitter.  Figure 
4(b)  shows  the  geometry  of  the  reduced  problem  where  the  scatterers  that  are  not 
close  to  the  receiver  and  transmitter  are  replaced  with  a  effective  dielectric  constant. 
Our  objective  here  is  to  show  that  the  mean  and  variance  of  the  field  in  the  original 
problem  and  the  reduced  problem  are  the  same.  The  method  of  moments  is  used  to 
solve  these  problems  for  a  given  arrangement  of  cylinders.  Then  using  Monte-Carlo 
simulation,  the  desired  field  statistics  are  computed.  In  the  implementation  of  the 
Monte-Carlo  simulation,  the  boundaries  of  the  medium  must  be  varied  randomly  by 
a  few  wavelengths  to  avoid  coherence  effects  [4],  Position  of  cylinders  are  determined 
by  a  random  number  generator.  In  this  filling  process  the  distance  of  a  new  cylinder 
is  measured  from  the  previous  ones  to  ensure  a  minimum  distance  between  the  tree 
trunks.  The  filling  process  for  each  medium  realization  is  continued  until  a  desired 
number  density  is  reached. 

3.1  The  Method  of  Moments  for  M-body  Sparse  Scatterers 

In  this  section  a  numerical  technique  based  on  the  method  of  moments  appropriate 
for  a  relatively  large  number  of  sparsely  distributed  dielectric  cylinders,  illuminated 
at  oblique  incidence  is  described.  The  integral  equation  formulation  for  the  induced 
polarization  current  for  infinite  cylinders  whose  axes  are  parallel  to  the  z-axis  is  given 
by 

1  N  r 

- - TJ (P)  =  ~ jhYoE i  +  £  /  G(\p-  pi)  •  J {(f)ds'  (6) 
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where  El  represents  the  incident  field  having  a  propagation  constant  k0,  er  is  the 
relative  dielectric  constant  of  the  cylinders  and  s'n  is  the  cross  section  of  the  nth 
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cylinder.  The  explicit  expression  for  the  dyadic  Green’s  function  is  given  by 


<5(1?-  ?\)  =  { 


k2  4-  —  &2  —  ik  — 

^0  ^  dx2  dxdy  JK'Zd x 

d2  J.2  .  SP_  _AU  d_ 

dxdy  ^0  1-  af  J"'z  dy 


H{J\kp\p-jj\) 


where  kz  is  the  propagation  constant  of  the  incident  wave  along  the  z-axis,  and  kp  = 
y/kQ  —  k2.  Following  the  standard  procedure  of  the  MoM,  the  cross  section  of  the 
cylinders  are  descretized  into  small  cells  over  which  the  current  distribution  may  be 
considered  a  constant  vector.  Using  point  matching,  the  integral  equation (6)  can  be 
cast  in  terms  of  a  matrix  equation  of  the  form 


Z- J  =  V 


(8) 


where  Z  is  known  as  the  impedance  matrix.  The  size  of  the  impedance  matrix  is 
proportional  to  the  total  number  of  cells  for  all  N  cylinders  and  is  a  limiting  factor 
with  regard  to  computer  memory.  Noting  that  in  a  fully-grown  forest,  the  tree  trunks 
are  sparsely  distributed,  storage  of  all  elements  of  Z  can  be  avoided.  In  this  case,  the 
impedance  matrix  of  the  individual  cylinders  (block  diagonal  elements)  are  computed 
and  stored.  This  can  be  limited  to  a  few  matrices  corresponding  a  small  number  of 
cylinder,  with  different  radii  that  represents  the  variability  in  tree  trunk  diameters. 
Next  the  impedance  matrix  elements  or  the  pairwise  interaction  between  the  cells 
at  the  center  of  the  cylinders  are  computed  and  stored.  The  interaction  between 
other  elements  are  computed  in  an  approximate  manner  as  needed  in  the  program 
and  are  not  stored.  For  further  clarification,  consider  the  ith  and  the  jth  cylinders 
in  the  global  coordinate  system  whose  centers  are  respectively  located  at  ( Xi,Yi )  and 
( Xj,Yj )  as  shown  in  Fig.  5.  Suppose  the  interaction  between  two  cells,  whose  local 
coordinate  in  the  ith  and  jth  system  are  respectively  given  by  (a;'///')  and  (a;'-,?/'),  are 
needed.  In  this  case,  using  the  Taylor  series  expansion, 


I  Pi  Pj  I  ~  Pij  T 


(Xi  —  Xj)(x'i  —  x'j)  _  (Xi-YjWt-it) 


Pij 


+ 


Pij 


(9) 


where  plj  =  y/(Xi  -  Xj)2  +  (T*  -  Y))2  is  the  distance  between  the  centers  of  the 
cylinders.  The  approximation  in  (9)  can  be  used  in  the  evaluation  of  integrals  needed 
for  the  computation  of  matrix  elements.  For  example, 


Ii'j'  — 


rx'j+A/2  />*/'+ A/2 

J  l' —  A/2  J  y'_  —  L 


H^ik.lpi  -  pj\)dx'dy' 

■i  -*j  A/2 

Wlu  ~  \  —  t 


anH{0  (kpPij)  ■ 
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where  §Vj>  =  kp  [(X*  -  Xj)(sJ  -  x'3)  +  (*i  -  Yj){y[  -  ?/'•)]  /  pi3,  and  an  =  A2  1  - 
As  mentioned  earlier  Ii3  is  stored  and  Ii>j'  is  computed  as  needed  .  Using  a  sparse 
matrix  storage  scheme  [5]  and  storing  a  small  number  of  terms  like,  (X*  —  X  j)/pij, 
(Yi  —  Yj)/pij ,  (a:'  —  x'j)  and  ( y[  —  ?/'•),  the  needed  memory  size  is  reduced.  Similarly 
for  other  terms  of  the  dyadic  Green’s  function  that  require  spatial  derivatives,  we  can 
use  the  following  approximations: 


du 2 


+  il  li’j1  ~ 


Oinkl 


(kpPij)  ±  H^\kppij)[ cos2  Oily  -  sin2  6??] }  -e^'  (10) 


where  the  positive  and  negative  signs  are  used  for  u  =  x  for  u  =  y,  respectively.  Also 
in  (10) 

COS2  evi,  -  Sin2  6Vf  »  cos2  ^-sin^^  +Cai-aal/pij 
3  13  l+(ai+a2)/Py  ’ 

at  =  2(X,  -  X,)(x'  -  f),  a2  =  2(y  -  Yj)(y '  -  y])  %  =  tan"1  [(y  -  Y^(X,  -  X)j 
Moreover 


dli’  <ynkl  r 


Ho]  (kppij)  +  H^ikppij)  I  (up  -  Ilf):  e*<v 


r(i) 


du  2 

for  u  —  x  or  y.  Finally  for  the  term,  we  can  use 

,,  =  (X,  -  x,)(y  -  y)  +  (X,  - x,)fo;  -  yj)  +  (y  -  y)(x-  - t_2 „(1)^  _ , 

dxdy  p?j  +  (qj  +  02) 


•kpH^ikppi^-e**'* 


Table  1  shows  the  sub-matrices  and  arrays  that  are  stored  in  this  scheme  for  M  iden¬ 
tical  cylinder,  each  descretized  into  N  pixels. 


3.2  Statistically  Equivalent  2-D  Medium 

In  this  section,  a  complete  MoM  solution  for  a  large  number  of  cylinders  is  used  to 
verify  the  conjectures  mentioned  for  the  reduced  problem.  However  before  presenting 
these  results,  the  accuracy  of  the  approximate  MoM  solution  for  the  M-body,  sparse 
scatterers  problem  is  evaluated.  An  array  of  50  equally  spaced  dielectric  cylinders 
with  er  —  5  +  j  and  radius  0.3m  are  arranged  along  Y-axis  with  a  spacing  of  2m. 
This  array  is  illuminated  by  a  50MHz  plane  wave  at  an  oblique  incident  angle  6  = 
60°  and  TE  polarization.  The  scattered  field  computed  by  the  exact  MoM  and  the 


approximate  MoM  are  compared  along  the  curve,  y  =  rrlnz.  Small  discrepancies  are 
observed  between  the  two  solutions  which  are  plotted  in  Fig.  7  as  relative  percentage 
error.  Many  other  cases  were  also  examined  and  it  was  found  that  when  the  average 
ratio  of  cylinder  spacing  to  cylinder  radius  is  larger  than  5,  the  approximate  MoM 
is  capable  of  producing  very  accurate  results.  Having  confidence  in  the  approximate 
MoM  algorithm,  simulation  of  wave  propagation  in  a  sparse  random  medium  is  carried 
out.  A  large  number  of  dielectric  cylinders  (2000)  confined  in  a  rectangular  box  as 
shown  in  Fig.  4,  are  considered.  The  box  is  chosen  to  have  average  dimensions  of 
50m  x  800m.  As  mentioned  earlier,  to  eliminate  the  coherence  effect  of  finite  box 
size  (see  [4])  the  dimension  of  the  box  is  varied  by  at  least  one  wavelength  in  the 
Monte  Carlo  simulation.  A  line  source  is  used  as  the  transmitter,  at  location  (x  = 
20,  y  =  25)  and  the  field  is  calculated  at  a  point  located  at  (x  =  780,  y  —  25). 
Many  sample  media  are  generated  and  the  MoM  solution  for  each  sample  is  used  to 
construct  approximate  statistics  of  the  field  at  the  receiving  point.  In  specifics,  the 
mean  and  standard  deviation  (SDV)  of  the  field  is  monitored  and  the  simulation  is 
continued  up  until  a  convergence  is  reached.  The  contribution  to  the  total  mean- 
field  and  standard  deviation  from  scatterers  in  range  bins  of  20m  wide  are  stored 
separately  and  are  plotted  in  Fig.  8,  for  three  frequencies  of  10,  30,  and  50MHz. 
It  is  interesting  to  note  that  the  contribution  to  the  total  standard  deviation  (field 
fluctuations)  are  mainly  dominated  by  the  scatterers  near  the  source  and  observation 
points.  This  indicates  that  scatterers  between,  but  not  close  to,  the  transmitter  and 
receiver  are  not  significant  sources  of  field  fluctuations  and  can  be  replaced  with  an 
effective  dielectric  constant.  The  geometry  of  the  reduced  problem  is  shown  in  Fig. 
4(b)  where  600  cylinders  near  the  source  and  250  cylinders  near  the  observation  point 
are  kept  and  the  rest  are  replaced  with  an  effect  dielectric  medium.  In  this  simulation 
a  cylinder  density  0.05/m2  is  used  and  the  comparison  of  standard  deviation  for  the 
complete  and  reduced  problems  for  the  TM  case  is  shown  in  Fig.  9.  As  can  be 
seen  in  Fig.  9  the  reduced  and  complete  problems  produce  effectively  the  same  field 
variations.  To  examine  the  field  variance  only,  far  less  scatterers  produce  the  same 
result.  Figures  10(a),  10(b), and  10(c)  show  respectively  the  field  variance  for  TE 
(10(a),  and  10(b))  and  TM  (10(c))  excitation  using  200  cylinders  near  the  source  and 
150  cylinders  near  the  receiver.  For  these  simulations  cylinder  density  0.01/m2  is 
used.  For  the  reduced  problem,  the  effective  medium  is  anisotropic  and  its  dielectric 
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constant  is  a  tensor  given  by  [4], 


Seff  =  £h  +  (£i  ~  £h)f  for  TM 

2 

£eff  =  £h  +  f(£i-  £h) - : -  forTE 

£{  +  £h 


4  Fast  Field  Calculation  Based  on  Reciprocity 

As  our  goal  in  this  investigation  is  the  computation  of  the  fields  of  a  dipole  in  a  forest 
environment,  computation  of  polarization  currents  in  tree  trunks  when  illuminated  by 
the  dipole  is  needed.  Even  with  the  simplification  mentioned  in  the  previous  section, 
a  3-D  scattering  problem  which  includes  more  than  200  cylinders  is  not  tractable 
numerically.  In  this  section,  we  demonstrate  a  procedure  where,  with  the  help  of 
reciprocity  theorem,  the  3-D  scattering  problem  is  first  reduced  to  an  equivalent  2- 
D  problem,  which  is  then  solved  by  the  method  of  moments.  To  demonstrate  this 
procedure,  let  us  first  consider  a  short  dipole  near  a  single  cylinder  as  shown  in  Fig. 
6.  The  field  at  the  receiver  is  the  sum  of  the  field  of  the  dipole  and  the  radiated 
field  from  the  polarization  current  induced  in  the  dielectric  cylinder  embedded  in 
the  canopy  above  the  ground.  Since  the  observation  point  is  in  the  far-field  of  the 
cylinder  and  the  dipole,  reciprocity  can  be  applied  to  simplify  the  problem.  According 
to  the  reciprocity  theorem  [6],  the  vertical  component  of  the  received  field  for  a 
dipole  excitation  with  orientation  l  is  equal  to  the  l  component  of  the  field  near  the 
cylinder  of  the  same  dipole  oriented  vertically  and  located  at  the  observation  point. 
According  to  (1),  the  field  of  the  dipole  (in  the  modified  problem)  illuminating  the 
dielectric  cylinder  is  locally  plane  wave.  Also  noting  that  the  induced  polarization 
currents  in  a  finite,  long  dielectric  cylinder  is  approximately  the  same  as  those  of 
an  infinite  cylinder  having  the  same  radius  and  dielectric  constant  [7,  8],  the  MoM 
solution  for  2-D  problems  can  be  used  to  find  the  induced  polarization  currents.  Once 
this  polarization  current  is  obtained,  the  near  field  can  easily  be  computed  and  the 
expression  for  it  given  in  the  Appendix.  The  same  procedure  is  applied  to  find  the 
contribution  of  all  cylinders  in  the  vicinity  of  the  transmitter,  at  the  receiver.  To 
compute  the  effect  of  the  scattered  field  of  cylinders  near  the  receiver,  the  fields  of 
the  dipole  and  all  cylinders  in  its  vicinity  are  computed  at  each  cylinder  location  near 
the  receiver.  Again  these  fields  are  locally  plane  waves  illuminating  tree  trunks  near 
the  receiver  at  an  oblique  incidence  equal  to  the  critical  angle.  MoM  is  used  to  find  the 
induced  polarization  current  in  cylinders  near  the  receiver  from  which  the  scattered 
field  is  computed.  Hence  the  contribution  from  all  cylinders  near  the  transmitter 


and  receiver  are  included.  To  account  for  the  effect  of  the  forest  floor,  the  geometric 
optics  images  of  the  dipole  and  lateral  waves  on  the  ground  plane  are  considered  and 
their  contributions  are  evaluated  using  a  similar  procedure.  It  is  worth  mentioning 
that  in  this  case,  the  number  of  cylinders  which  need  to  be  kept  in  the  vicinity  of  the 
transmitter  and  receiver  for  the  reduced  problem  is  expected  to  be  smaller  than  what 
was  obtained  for  the  2-D  problem.  Basically  for  the  dipole  excitation  where  the  finite 
tree  trunks  are  illuminated  by  the  resulting  plane  waves  at  oblique  incidence,  fewer 
number  of  cylinders  can  interact  with  each  other. 


5  Numerical  Simulation 

Based  on  the  rigorous  electromagnetic  model  described  in  the  previous  section,  a 
very  accurate  propagation  model  which  provides  a  complete  channel  characterization 
of  forest  media  is  possible  through  Monte  Carlo  simulations.  To  examine  the  effect 
of  tree  trunks  on  the  field  of  a  transmitter  in  a  forest,  we  first  consider  a  single  tree 
trunk  in  the  near-field  of  a  short  dipole.  For  this  example  a  dielectric  cylinder  with 
permittivity  e  =  5  +  j,  radius  a  =  35cm  ,  and  height  hc  =  15m  is  considered  in  a 
forest  with  a  canopy  having  £ejf  =  1.03  +  jr'0.036  and  height  H  =  20m.  A  vertical 
dipole  transmitting  at  90MHz  is  placed  at  a  distance  of  10cm  from  the  surface  of  the 
cylinder  and  is  moved  up  and  down,  and  around  the  cylinder  and  its  field  is  com¬ 
puted  at  an  observation  point  1km  away  from  the  cylinder  and  5m  above  the  ground 
plane.  Figure  11  shows  the  normalized  z-component  of  the  field  at  the  receiver  as 
a  function  of  dipole  height  above  the  ground  and  for  three  azimuthal  angles  around 
the  cylinder.  The  normalization  here  is  with  respect  to  the  field  of  the  dipole  in  the 
absence  of  the  cylinder.  It  is  shown  that  the  field  at  the  receiving  point  fluctuate  as 
the  dipole  is  moved  along  the  cylinder  axis.  This  fluctuation  is  a  result  of  construc¬ 
tive  and  destructive  interference  of  the  direct  field,  scattered  field  from  the  cylinder, 
and  their  images  on  the  ground  plane.  Since  the  cylinder  radius  is  small  compared 
to  the  wavelength  a  gentle  variation  is  observed  as  the  dipole  is  moved  around  the 
cylinder.  It  is  interesting  to  note  that  for  most  dipole  locations,  existence  of  the 
cylinder  in  the  near-field  region  of  the  dipole  enhances  the  field  at  the  receiver  since 
the  tree  trunk  acts  as  a  passive  radiator.  Next,  computation  of  path-loss  and  field 
standard  deviation  is  considered  for  the  forest  of  the  previous  example.  In  this  case 
tree  trunks  having  an  average  height  15m  and  height  standard  deviation  of  lm  with 
number  density  0.05/m2  and  dielectric  constant  of  e  =  5  +  j  are  also  included.  As 
mentioned  before,  the  number  of  cylinders  we  need  to  keep  for  the  reduced  problem 
is  expected  be  smaller  than  those  for  the  2-D  problem.  Here  we  kept  200  cylinders 
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near  the  transmitter  located  at  the  origin  3m  above  the  ground  and  50  cylinders  near 
the  receiver  located  1km  away  from  the  transmitter  and  5m  above  the  ground.  Table 
2  shows  the  path-loss  (with  respect  to  the  free-space)  of  the  normal  component  of  the 
wave  and  its  standard  deviation  at  50MHz  for  a  vertical  dipole.  Also  shown  in  Table 
2  are  the  path-loss  and  standard  deviation  if  390  cylinders  (250  near  the  source  and 
140  near  the  receiver)  are  kept.  The  Monte-Carlo  simulation  converges  after  about 
30  realizations  and  is  shown  that  the  results  based  on  200  cylinders  is  very  close  to 
those  based  on  390  cylinders.  This  convergence  test  indicates  that  a  moderate  number 
of  scatterers  (about  200)  is  sufficient  for  accurate  characterization  of  field  variance. 
The  results  also  indicate  a  standard-to-mean  ratio  deviation  of  almost  unity  for  the 
random  variable  Ez.  Using  different  forest  parameters  and  frequencies  (at  HF  band) 
random  variable  Ez  showed  a  similar  behavior.  That  is  the  statistics  of  the  channel 
follows  Ricean  statistics  with  a  standard  deviation-to-mean  ratio  ranging  from  0.5  - 
2.  This  is  different  than  Rayleigh  statistics  which  is  commonly  assumed  for  commu¬ 
nication  channels.  Significant  non-zero  mean-field  is  a  direct  result  of  contributions 
from  the  lateral  wave. 

For  communication  channels  with  significant  multi-path,  antenna  arrays  are  usually 
used  to  mitigate  fading.  In  this  case,  because  of  the  existence  of  a  considerable  co¬ 
herent  mean-field,  antenna  arrays  may  provide  coherent  gain.  To  investigate  the  per¬ 
formance  of  antenna  arrays  in  a  forest  environment  two  basic  configurations,  namely, 
broadside  and  end-fire  are  considered.  For  the  broadside  configuration  four  vertical 
dipoles  are  arranged  along  a  line  perpendicular  to  the  line  between  the  transmitter 
and  receiver  points  with  a  spacing  of  A/2.  The  four  elements  of  the  end-fire  array 
are  placed  along  the  line  between  the  transmitter  and  receiver  separated  by  7A/16 
and  having  a  progressive  phase  factor  of  -77r/8  [9].  The  field  of  these  two  arrays 
are  evaluated  and  the  path-loss  and  SDV-to-mean  ratio  is  reported  in  Table  2.  It  is 
found  that  the  broadside  array  has  about  lldB  less  path-loss  than  the  end-fire  array. 
This  is  very  close  to  the  12-dB  gain  of  a  4-element  array  and  indicates  that  the  field 
of  all  four  elements  of  the  broadside  array  arrive  at  the  receivers  almost  coherently 
whereas  those  of  the  end-fire  array  are  incoherent.  This  behavior  can  be  understood 
by  considering  the  spatial  correlation  function  in  the  forest  environment.  The  spatial 
correlation  coefficient  between  two  points  fi  and  r*2  is  defined  as  [11] 

Cov  \E(fl),E(r2) 

Cs(r1,f2 )  = - - 

where  oei  and  <je2  are  the  variance  of  the  field  at  r)  and  r2  respectively.  A  grid  of 
nine  points  along  x-axis  (direction  between  the  transmitter  and  receiver)  and  nine 


points  along  y-axis  separated  by  A/2  that  is  1km  from  the  transmitter  is  generated 
and  components  of  the  electric  field  are  calculated  many  times  for  different  cylinder 
arrangements  in  order  to  compute  the  correlation  coefficient.  Figures  12(a)  and  12(b) 
show  the  magnitude  and  phase  of  Cs(fi,  rjj)  along  the  x-  and  y-axis  respectively.  It  is 
shown  that  the  signal  along  y-axis  remains  coherent  over  much  larger  distance  than 
along  x-axis.  This  result  is  in  agreement  with  the  observed  behavior  of  broadside  and 
end-fire  antenna  arrays.  Sampling  theorem  [10,  11,  12]  is  used  to  generate  the  smooth 
curves  that  goes  in  between  the  discrete  points  in  Figs.  12. 

To  study  the  path-delay  effects,  the  impulse  response  of  the  channel  must  be  char¬ 
acterized.  The  impulse  response  cannot  be  directly  determined.  However,  by  char¬ 
acterizing  frequency  response  over  a  wide  bandwidth  the  impulse  response  can  be 
determined  by  applying  the  Fourier  transformation  to  generate  time  domain  (range) 
information.  Direct  application  of  FFT  is  not  computationally  efficient  considering 
the  number  of  frequency  points  which  are  needed  in  order  to  avoid  aliasing.  For 
example,  with  a  receiver  at  a  distance  of  1.5km  from  the  transmitter  a  maximum 
frequency  spacing  of  200  KHz  is  required.  A  pulse  width  of  12.5  nsec  corresponds  to 
80MHz  bandwidth  which  can  be  used  to  resolve  path  delays  3.75m  apart.  In  this  case 
for  each  forest  realization  400  frequency  points  must  be  simulated.  To  circumvent 
this  difficulty  we  used  a  non-uniform  frequency  sampling  scheme.  Gauss  quadrature 
integration  [12]  is  used  to  evaluate  the  Fourier  transform.  The  order  of  Legendre 
polynomial  should  be  chosen  so  that  the  minimum  distance  between  the  zeros  of  the 
Legendre  polynomial  is  smaller  than  the  minimum  frequency  spacing  required  to  avoid 
aliasing  (200  KHz  for  the  above  example).  Impulse  response  of  the  forest  considered 
in  the  previous  example  is  simulated  at  HF  using  a  bandwidth  of  10  -  90MHz.  A 
receiver  at  a  distance  of  1km  from  the  transmitter  is  considered.  It  is  expected  that 
the  dispersion  (pulse  spreading)  will  be  observed  due  to  multiple  scattering  among 
tree  trunks  and  the  frequency  dependent  path-loss  of  lateral  waves.  Figure  13  shows 
the  frequency  response  of  the  received  field  in  the  absence  of  the  tree  trunks  and  the 
ground  plane  for  two  cases:  1)  smooth  canopy-air  interface  and  2)  rough  canopy-air 
interface  with  a  rms  height  2m.  The  transmitter  and  receiver  are  17m  and  15  m  below 
the  interface  and  the  effective  dielectric  constant  of  the  canopy  is  eefj  —  1. 03 +j  0.036, 
as  before.  It  is  shown  that  the  canopy-air  interface  roughness  increase  the  dispersion. 
Figure  14  shows  the  impulse  response  of  the  forest  channel  including  tree  trunks  and 
the  underlying  ground  plane.  For  the  calculation  of  the  impulse  response,  Gauss 
quadrature  method  with  40  points  is  used  over  the  frequency  range  of  10  -  90MHz. 
A  Gaussian  pulse  is  assumed  to  be  transmitted  which  is  also  plotted  in  Fig.  14  for 
comparison.  The  amplitude  of  the  transmitted  pulse  is  multiplied  by  the  path-loss 
and  delayed  by  the  free-space  distance  between  the  transmitter  and  receiver.  As  can 
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be  seen  in  Figure  14,  pulse  spreading  and  ringing  are  observed  which  are  the  result  of 
dispersion  and  multipath.  The  last  simulation  demonstrates  performance  of  different 
spatial  diversity  schemes.  Three  antennae  1.5A  apart  are  aligned  x-axis  placed  1km 
away  from  a  transmitter  operating  at  50MHz  in  the  forest.  Three  diversity  schemes 
are  examined:  1)  selective  diversity  (SD),  2)  equal  gain  combining  (EG),  and  3)  maxi¬ 
mal  ratio  (MR)  combining  [13].  In  SD  scheme  the  detector  simply  chooses  the  output 
of  the  receiver  with  the  highest  receiver  power.  In  EG  combining  method  the  detector 
is  provided  with  the  averaged  detected  signal  from  each  receiver  (a  weighting  factor  of 
Wi  =  1/3  is  used  in  the  combining).  In  the  diversity  scheme  based  on  maximal  ratio 
combining  a  weighting  factor  proportional  to  the  signal  power  is  used  to  combine  the 
signals  from  the  receiver  (wi  —  ).  To  assess  the  performance  of  the  above- 

mentioned  diversity  combining  methods,  cumulative  distribution  function(CDF)  of 
fading  depth  for  each  diversity  scheme  is  calculated  using  a  Monte-Carlo  simulations. 
The  fading  depth  is  defined  as  power  level  at  the  output  of  each  combiner  (SD,  EG, 
or  MR)  above  the  average  power  in  dB  (101og(|i?|2/(|.E|2))).  Figure  15  shows  the 
CDF  of  fading  depth  for  SD,  EG,  and  MR  combining  methods.  Also  shown  is  the 
CDF  of  one  of  the  channels.  It  is  shown  that  the  performance  of  all  three  diversity 
scheme  is  approximately  the  same  for  this  problem. 

6  Concluding  Remarks 

A  complete  wave  propagation  model  capable  of  predicting  path-loss,  time  delay  re¬ 
sponse,  and  coherent  frequency  response  in  a  forest  environment  is  developed.  The 
model  is  based  on  a  hybrid  analytical  and  numerical  approach  which  allows  for  effi¬ 
cient  computation  of  important  channel  properties  without  compromising  accuracy. 
The  model  accounts  for  multiple  scattering  among  tree  trunks  and  includes  the  effects 
of  the  forest  ground  plane.  Branches  and  leaves  are  represented  by  an  effective  dielec¬ 
tric  constant  which  limits  the  range  of  validity  of  the  model  up  to  UHF.  Except  for 
the  path-loss  calculation  other  statistical  channel  characteristics  are  determined  using 
a  Monte-Carlo  simulation.  In  general,  propagation  simulations  are  slower  at  higher 
frequencies  which  depending  on  the  computer  platform  also  constitutes  a  limit  for 
the  highest  frequency.  The  model  is  used  to  simulate  performance  of  antenna  arrays 
and  different  diversity  schemes.  It  is  shown  that  because  of  a  significant  non-zero 
mean-field,  spatial  diversity  only  improves  the  channel  fading  moderately  and  there 
is  not  a  significant  difference  between  different  diversity  combining  schemes.  The 
model  presented  in  this  paper  shows  the  possibility  of  constructing  a  very  accurate 
physics-based  propagation  model  capable  of  end-to-end  channel  simulations. 
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Appendix:  MoM  Formulation 

In  this  appendix  a  general  method  of  moments  formulation  for  2-D  dielectric  objects 
with  possible  dielectric  inhomogeneity  is  provided.  The  axis  of  a  cylinder  with  an 
arbitrary  cross  section  is  assumed  to  be  parallel  to  the  z  axis  and  the  surrounding 
medium  is  assumed  to  be  free-space.  Let  the  relative  permittivity  of  the  cylinder  be 
er(x,y)  and  its  relative  permeability  be  unity  (pr  =  1).  The  cylinder  perturbs  the 
incident  field  and  the  difference  between  the  perturbed  (total)  and  incident  field  is 
known  as  the  scattered  field;  thus, 


E‘(i>)  =  E‘G<>)  +  E'(/>), 


(A.l) 


where  p  is  the  position  vector  in  cylindrical  coordinates.  From  Maxwell’s  equations 
it  can  be  shown  that  a  volumetric  current  density  of  the  form 


J e(p)  =  - ik0Y0[er(p )  -  1]E4 (p), 


peS, 


(A.2) 


known  as  the  polarization  current,  in  free  space,  can  replace  the  cylinder  to  reproduce 
the  scattered  field.  Therefore  the  scattered  field,  can  be  obtained  from: 


\p)  =  -ik0Z0  [ Je(p')-G(p,p')ds, 

J  s 


(A.3) 


where  G(p,p)  is  the  two-dimensional  dyadic  Green’s  function  given  by  ([7]).  Using 
(A.l),  (A. 2), and  (A.3)  an  integral  equation  for  the  unknown  polarization  current  can 
be  obtained, 

J e(p)  =  -ik0Y0[er(p)  -  lKE^p)  -  ik0Z0J^e(p')-^(p,p,)ds}. 

Assuming  the  incident  field  is  illuminating  the  cylinder  at  an  oblique  incident  angle 
Oi  (k\  =  k0  cos  Qi),  in  order  to  satisfy  the  phase  matching  condition,  all  field  and 
polarization  current  components  have  a  z-dependance  of  the  form  e~kzZ.  If  the  incident 
field  is  normal  incident  klz  =  0,  then  the  problem  may  be  decoupled  into  TM  and  TE 
problems.  For  the  TM  case  the  incident,  scattered,  and  hence,  the  polarization  current 
have  only  z  components  and  the  integral  equation  is  simplified  to 


Jz(x,  y)  =  —ik0Y0[er  (x,  y)  -  l]££(x,  y)  +  i^[er(x,  y)  -  1] 
fs  Jz(x',y')H^(kp  I  p-  p  |)dx'dy'}. 


(A.4) 
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In  the  TE  case  both  x  and  y  components  of  the  polarization  current  are  induced  and 
they  satisfy  the  following  coupled  integro-differential  equations 


Jx{z,y)  =  -ik0Y0[er(x,y)  -  l]J5£(a:,y)  +  f[er( x,y)  -  1] 

■{(&  +  ko )  Is  y')Ho1](kP  [p-  p  \)dxdy 

+&  Is  Jy(x'>  y')Ho](kP  I  p-  p  \)dxdy'}  .  > 

Jy(x,y)=  -ikoYo[^r(x,y)-l]E^x,y)  +  l\er(x,y)-l]  {  } 

■{-&iSsJx(x'’y')HQ  \kP  I  P-P 

+(£2  +  k%)  f$  Jy{x',  y')H(I\kp  Ip-  p  |)dx'dy'} 

The  resultant  integro-differential  operators  obtained  for  the  polarization  current  do 
not  impose  any  restriction  on  the  functional  form  of  the  current,  in  particular,  the 
pulse  function  is  in  the  domain  of  the  operators.  It  should  be  noted  here  that  the 
kernel  of  the  integral  equation  (Green’s  function)  for  the  TE  case  is  more  singular  (4) 
than  for  the  TM  case  (In  p) .  There  is  no  known  solution  for  these  integral  equations 
in  general,  but  their  forms  are  amenable  to  numerical  solution.  An  approximate 
numerical  solution  for  the  integral  equations  can  be  obtained  using  the  method  of 
moment  in  conjunction  with  the  pulse  expansion  function  and  the  point-matching 
technique.  The  cross  section  of  the  scatterer  is  divided  into  N  rectangular  cells  that 
are  small  enough  so  that  the  polarization  current  and  relative  permittivity  can  be 
assumed  to  be  constant.  The  unknown  current  can  be  approximated  by 

N 

Jp{x,y)  =  ^2jmP(x-xm,y-ym),  p  =  x,y,orz  (A.6) 


where  Jm  are  the  unknown  coefficients  to  be  determined  and  P(x  —  xm,y  —  ym)  is  the 
pulse  function  defined  by 


P{x  xm,  y  ym) 


1  |  X-Xm  |<  \V-Vm\<^$ 

0  otherwise 


(A.7) 


By  inserting  the  current  as  expanded  in  (A.6)  into  the  integral  equations  (A.4)  and 
(A. 5)  and  then  setting  the  observation  point  at  the  center  of  the  mth  cell,  a  linear  set 
of  equations  is  formed.  In  matrix  notation,  these  linear  equations  can  be  represented 
by 

Z™Jz  =  £z  (A. 8) 

for  the  TM  case,  where  Z™  is  the  impedance  matrix,  Jz  is  the  unknown  vector,  and 
[£z]  is  the  excitation  vector.  Similarly  for  the  TE  case  the  coupled  integral  equation 
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(A. 5)  in  matrix  form  becomes 


ZEEJX  +  ZEEJy  =  £x 

■p'TE  rr  ,  jrTE  q  _  c 
"3  Ox  I  "4  Jy  —  C-y 

where  as  before  ZfE,  ■  ■  ■  ,  ZjE  are  TV  x  TV  impedance  matrices  and  Sx  and  £y  are 
the  excitation  vectors.  The  above  coupled  matrix  equation  can  be  represented  by  a 
27V  x  27V  matrix  equation  similar  to  (A. 8).  For  the  general  case  of  obliquely  incident, 
all  three  current  components  are  coupled,  and  the  final  matrix  becomes 


ZEEJX  +  ZEEJy  +  ZxJz 

ZEEJX  +  zJEjy  +  zyjz 

Zx.7„  4-  Zy.7..  +  k2Z™. 7. 


(A.10) 


Although  the  variation  of  the  polarization  current  and  dielectric  constant  over 
each  cell  is  ignored,  this  cannot  be  done  for  the  Green’s  function.  Actually  for  cells 
close  to  the  observation  point,  the  Green’s  function  varies  rapidly  and  its  contribution 
must  be  evaluated  more  precisely.  Let  us  denote  the  function  representing  the  Green’s 
function  contribution  by 

fXn  +  ±t  rVn  +  ^  _ 

In(x,y)=  /  Hi  \kpyj{x  -  x'Y  +  {y-  y')2)dx'dy'  (A.ll) 

Jx  »-£f*  Jyn-^ 


where  1  <  n  <  TV.  If  the  observation  point  (x,  y)  is  different  from  ( xn,yn )  ,  then  the 
integrand  in  (A.ll)  is  not  singular  and  since  )  x'  —  xn  |<  A2G.  an(j  |  y'  —  yn  |<  its 
Taylor  series  expansion  may  be  substituted.  By  retaining  the  terms  up  to  the  cubic 
order  in  the  expansion  of  the  Hankel  function,  In(x,y)  is  found  to  be: 


In(x,  y)  =  AXnAYn{H^] {kpy/(x  -  xn)2  +  {y-  yn )2)  +  (fc^fn)2  A(x  -xn,y-  yn)+ 
(-p*Jn)2 B(x  -xn,y-  yn)}, 

(A.12) 

where 


A(x  -  xn,y-  yn)  =  A(rn,  9n )  =  -H§\kprn)  cos2  9n  + 


B(x  -  xn,y-  yn)  =  B(rn,  9n)  =  -H^\kprn)  sin2  9n  + 


Hx  (kprn)  (cQg2  _  gin2 

V" 

(A.13) 

Hi  (kprn)  ^.n2  _  cqs2 

(A.14) 
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with  the  following  definition  for  a  pair  of  local  polar  coordinates: 


rn  =  ^{x  -  xn)2  +  (y  -  yn)2 
en=  arctan(l^). 


(A. 15) 


The  first  order  derivatives  of  In(x,  y)  are  also  easily  calculated  and  are  given  by 


dln(x,y ) 

fix 


=  AXnA Yn{~H^\kpy/{x  -  xny  +  {y-  yn)'2 )  cos  9n+  (A _16) 

-  xn,y  -  yn)  +  ~  *»>  y  ~  »»)> 


=  AXnAYn{-H[1]  ( kp  s/(x  -  xnf  +  (y  -  yn )2)  sin  9n+ 


[X  -  xnj~  -r  yy  -  y  am /A  17) 


where 


^(r„A)  =  H[l\k„ r„)cosA  -  ^1>(V»)as£^r!li‘+ 

/i V  rr(^) 

^^(cos2  -  sin2  0n)  A:pHi  (&prn) - J- 


g^fcpTn) 


H^jkpTn)  _  4  cos  fin,  sin2 

kpTn  Tn 

>,,»,)=  ifWtVOcos3^  - 

^(cosA-sinA)  [fcpffW'tV.)  -  StVii  - 

H^jkpTn)  _  4 cos  fln  sin2  Bn 
kprn  T n 

%A(rn,6n)  =  H?\kprn)  cos2  0n  sin  $n  +  H^(kprn)2cos2 

^^(cos2  en  -  sin2  en)  kpH x(1)  ( kprn )  -  AlJAsIA  + 

H[l\kprn)  _  4  cos2  sin  6n 
kpTn  Tn 

iB(r„#,)  =  //"' cos2  sin  ft,  + 

|^(cosA-sin20„)  kX‘VM  ~  ~ 

H^jkpTn)  _  4 cos2  6n  sin  9m 

kpTn  Tn 


(A.18) 


(A.19) 


(A. 20) 


(A.21) 
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The  second  order  derivatives  of  In(x,  y )  are  also  needed  for  calculation  of  the  impedance 
matrix  elements  for  the  TE  case,  and  are  given  by 


(A  +  k2)In(x,  v)  =  ’(Vw)  sin2  K  +  6,  -  sin2  »„) 

+^(&  +  «.)  +  +  *«2)B('-„,e„)}+ 

AX„A  Ynk]H^(kpr,)  sin2  9n 

(&  +  k2)Ux,y)  =  AX„AY„k2p{H^\kpTn)  cos2  fl„  +  (sin2  0n  -  cos2  #„) 


+  k20)A(rn,en)  +  ^(^  +  *S)B(r„,ft,)}+ 
AX„AY„k2H^(k„rn)  cos2  9n 


(A. 22) 


where 

a2 


^A(rn,0n) 


^A(rn,  0n)  =  ^{#<0(  Vn)[cos2  9n(§  cos2  0n  +  |  sin2  0n) 

+  ($nT^(3  C0s2  ^  “  Sin2  ^ 

+^1)(V^)[^cos2^sin2^  -  |f^(3cos20n  -  sin2  0n)] 

+F^1)(V»)[-|  C0S4  +  0rp(-9COs2  6 '»  +  ^  ^ 

+tff)(A;,,rn)[|  cos2  0n(cos2  0„  -  sin2  0„)]}, 

(A. 23) 

0„(!  COS2  0n  +  |  sin2  On) 

+  \l°pfn)2  (CQs2  0n  -  3  sin2  0„)] 

cos2  *» sin2  *»  +  ifSjH3 sin2  -  cos2  6n)] 

+Hi1)(kprn)[-\ sin2  9n  cos2  9n  +  g^r(9  sin2  0„  -  cos2  9n)} 
+H^\kprn)[l  sin2  0„( cos2  0n  -  sin2  0„)]},  ^  ^ 

£rB(r„,  9n)  =  k2p{H^  (kprn)[  cos2  0B(§  sin2  0n  +  §  cos2  0B) 

+  20;-y£- (sin2  6>n  “  3cos26>„)] 

+i/F(Vn)[^(-4cos2  0n  +  sin2  0B)  +  f0^(3cos2  9n  -  sin2  0B)] 
+^1} (kprn)[-\  sin2  9n  cos2  9n  +  (9  cos2  0„  -  sin2  9n)\ 

+H(4l\kprn)[l  cos2  0„(sin2  9n  -  cos2  0B)]}, 


(A. 25) 


+H[1}  (kprn)[^r  sin2  9n  cos2  9n  -  4f°fJ3n  (3  sin2  6n  -  cos2  6n)] 


j^B(rn,  en)  =  k2p{H^ (kprn)[sm2  0n( §  sin2  9n  +  f  cos2  6n) 

+  2{CrnY  (3sin2  °n  -  cos2  9n )] 

rF 

+^1}(A;prn)[-|  sin4  9n  -  g^(9sin2  9n  -  cos2  9n)\ 
+H{41)(kprn)[l  sin2  6>n(sin2  0n  -  cos2  (9n)]}. 

a 2 


(A. 26) 

We  also  note  that  an  exact  analytical  expression  for  ~^In(x,y)  can  be  obtained 
without  using  the  Taylor  expansion  and  is  given  by 

d2-In(x,y)=  H{Q]{kpJ  (x-xn-  ^)2  +  {y  -  yj 


dxdy 


H{0l\kp^(x-xn-^)2  +  (y-yn  +  ^)2)- 
H(01] ( kp ^(x-xn  +  ^)2  +  (y-yn-  ^)2) + 
HfHkpy/  (x-xn  +  ^)2  +  {y-yn  +  ~[L)2) 


(A. 27) 


When  the  observation  point  is  in  the  center  of  the  cell  itself,  the  Taylor  series  expan¬ 
sion  cannot  be  used.  In  this  case  we  can  employ  the  small  argument  expansion  of  the 
Hankel  function,  i.e. 

r(i)/  s  -  2 x 


Hq  {x)  ~  1  -I  (In  —  +  7) 

7 r  z 

Then  at  the  center  of  the  cell  (self-cell  contribution),  we  have 

_  i4fk2pAXnA  Ynj_  in+3  ,  ,„,kp^/(AXnr+(&Yn)2 


n(xn,yn)=  f{^J^[7-^+ln(' 
,  (k£AXls2  arH.on/ ,45k\  .  (I 


)] 


+(^)2arctan(££)  +  (^)2(|  -  arctan(^))}. 


(A. 28) 


(A. 29) 


Using  the  same  expansion  we  can  also  get 


(£  +  kl)In(xn,yn)  =  7  -  M  +  ln(V^x„)^(A'-„) 


■)] 


...  [7  2 

+(^)2arctan(^)  +  (^)2(|  _  arctan(^))}, 


(A. 30) 


+  [7  -  +  in(^-AXf +(Ar")2)] 


+(^)2(f  -  arctan(^))  +  (^)2  arctan(^)}. 

(A.31) 

The  evaluation  of  the  second  order  derivatives  of  In(x,y )  (expressions  in  (A. 22)) 
gives  accurate  results  when  rn  >  A/60.  For  smaller  values  of  rn  the  small  argument 
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(A. 32) 


expression  for  the  Hankel  function  can  be  used.  In  such  cases  we  have 


(£2  +  kl)In(x,y)  =  Fi(x,  y)  —  F2(x,  y) 
(■§^  +  ko )In(x,y)=  Gi(x,y)  -  G2(x,y) 


where 


V )  =  -¥-!)  +  (‘an  gf  -  to"  tt )(“  - 

-  ^  Kj  In  i^S+5  -  On!  1„  ^^1 

*o24M“  -  ¥  -  1)  +  (<*”  !“  -  ‘an  £)(?  - 


(A. 33) 


IT 

io-nj  k2p  r u  1  y/^n2+a\ 

7 r  LUn2  111  2 


with 


i- 


+ 


■6n  iln^2+<] 

Ma  i  -  1 

aL  --2 

(A. 34) 

ASk  j  =  1 

4  2  =  2 

(A.35) 

Un 

Vn  + 

Now  we  are  in  a  position  to  express  the  impedance  matrix  elements  in  terms  of 
In(x,y).  The  off-diagonal  entries  of  the  impedance  matrix  for  the  TM  case  are  given 
by 

ikl. 


7tm 


_  [6y  (Xff 2)  Vm)  1  \^n{xm^ym) 


and  the  diagonal  entries  are 


Zln  =  *-~[er(xn,yn)  ~  1  ]/«(*«,  lft») 


(A. 36) 


(A. 37) 


For  TE  polarization,  where  the  impedance  matrix  is  composed  of  four  sub-impedance 
matrices,  the  off-diagonal  elements  of  each  matrix  are 


yTE  _ 
^lmn 
yTE  _ 

yTE  _ 

yTE  = 
^ 4mn 


4  [G  {%m  ■> 
4  [G  (^mi 

yTE 

“2mn 

4 


2/m)  ^][(aj2  T  k'o)In{xrn,  ym)] 

Vm)  ~  ^\\.QxQy^n{Fmi  2/ro)] 

2/m)  —  1] [( A  ko)In(xm ,  2/m)] 


and  the  diagonal  elements  are  given  by 

=  |MZn,yn)  “  l][(J^  +  kl)In{xn,  yn)]  -  1 

=  0 

Zjrm  ~  0 

ZInn  =  I [tr(Xn,  Vn)  ~  +  kl)In(xn,  yn )]  -  1 


(A. 38) 


(A. 39) 


Finally  the  terms,  Zx,  and  Zy  are  given  by 


yx 
^ mn 


kz 

4 


—mn  4  [er  2/ra)  l]  Vm) 

Z mn  =  f  Vm)  ~  1  ]g^kn(xm,ym) 


for  off-diagonal  elements,  and 


(A. 40) 


Zn n  =  Zynn  =  0  m^n  (A.41) 

for  diagonal  elements.  The  excitation  vector  elements  for  the  TM  and  TE  incidence 
cases,  respectively,  are  given  by 


[er(2;m,ym)  l\Ex(xm,  ym)i 

TE 

(A. 42) 

(^m>  Um)  ^-\Zy{xm,  ym) j 

TE 

(A.43) 

(xm,  ym)  l]£/z(xm,  ym) 

TM 

(A. 44) 

Appendix  B:  Near  Field  Calculation 


In  this  section  an  efficient  formulation  for  the  calculation  of  scattered  fields  from 

dielectric  cylinders  is  provided.  With  the  help  of  the  Hertz  potential,  electric  field 

— »  —* 

can  be  easily  computed  from  a  known  current  distribution,  J,  by  using  E  =  V(V  • 
ne)  +  k\YLe.  The  electrical  Hertz  vector  is  represented  by 


4?r  jwe  Jv, 


_  ...  ,e3^-r'\ 

Jne~]k‘z '  ^.dv' 


(B.l) 


where  Jn  is  a  constant  vector.  Let  /  =  fs,  1  e  ik*z'dx'dy'  =  fs,  /(|r— r'De  ik'zZ> dx' dy' 

.  where  s'n  is  a  cell  area.  For  small  pixel  area  and  using  mid-point  integration 


22 


where  rn  =  y/(x  —  xn)2  +  (y  -  yn )2  +  (z  -  z')2.  Based  on  the  above  result,  other 
integrals  can  be  evaluated  as 


I xx  =  a-?/  ~  A2 - (jfcrn  -  1)  — 


1  3  cos2  9n 


+  jk  cos2  #r 


r  5  r  2  ejAT"  . .  1  3  sin2  #n  .  2. 

Jw  =  - Ofcr„-l) - +  ^sm20n 

dy2  rn  rn  rn  J 


£*1  «  A- - (jkrn  -  1) 

<9z2  rn 


'  1 


3  COS2  9n 


+  jk  cos2  9r 


T  -  9  T~  A2jkrn  -  lr^  f  COS0n  M  = 

u  -  du  r2  \  sin0„  u  = 


The  rest  terms,  Ixy,  Ixz  and  Iyz,  can  be  evaluated  analytically  as 


Ixy  =  A2 


eJkrn  ejkrn+  ejkr+  ejkr++ 


cP-  rL  pj kr no  neikrnl  •  1 1  r 

Ixzdz'  =  — —  /  Idz'  =  A2 - ( jkrn0  -  1)  cos  9n0  -  A2 - ( jkrnl  -  1)  cos  9rde  jk* 

dxdz  J o  rn0  rnn 

Fp1  pjkrnQ  l/r 

Iyzdz'  =  — -  /  Idz'  =  A2 - (jkrn0  -  1)  sin0nO  -  A2 - (jkrni  -  1)  sm9nie~3  - 

dydz  J0  rn0  rnl 


where 

rn±  =  V(x~xn  ±  A/2)2  +  (y-yn±  A/2)2  +  (z  -  z')2 
rn 0  =  V(x~  xn)2  +  (y-  Vn)2  +  Z2 
rnl  =  V(x~  xn)2  +  (y  -  Vn)2  +  (z~  L)2 
c°S0„o  =  Z/rn0,  COS0nt  =  (z  -  L)/rn0 

Therefore  final  expression  for  the  electrical  field  that  includes  the  effect  of  ground  is 
given  by 

[  { (Rh  +  1  )(IXX  “b  k1I)lx  T  (Rh  T  P)Ixyly  T  (Ry  T  1  )Ixz^z} 

{(Rh  +  P)IXylX  +  ( Rh  +  ^){dyy  +  k\l)ly  +  (Ry  +  l)Iyzlz}  y+ 

{(Rh  +  1  )lxzlx  +  (Rh  +  1  )Iyzly  +  (Rv  +  1  )(^zz  +  ^1-0^}  + 

RhihzX  +  IyzV  +Izzz)( cos  (plx  +  sin  (ply)]  dz' 

(B-2) 
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where  the  summation  is  over  all  cells  in  all  cylinders  and  Rh,  and  Ry  are  the  ground 
Fresnel  reflection  coefficient  and  also  Rh  =  sin  29  cos 6  [14]. 

k  cos  6+ y  ft— sin2  6 
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Stored  Components 

Size 

Imp.  Matrix  of  individual  cylinder 

M  xM:  TM,  2 M  x  2 M:  TE 

kp(\i  X j ) / Pi j ,  kp(\ i  Yj)/pij 

M(M  -  1)  array 

$i'j' 

M(M  -  1)  array 

cos2  6ij  —  sin2  6ij 

M(M  -  l)/2  array 

Table  1:  Stored  matrices  and  arrays  for  the  implementation  of  M-body  sparse  scat¬ 
tered. 


Antenna  Configuration 

Path-loss  (dB) 

SDV /Mean(dB) 

200  cyl. 

390cyl. 

200  cyl. 

390  cyl. 

single  dipole 

47.5 

48.16 

1.42 

1.16 

4-element  Broadside 

48 

3 

4-element  End-fire 

59 

2.45 

Table  2:  Table  2:  Path-loss  and  standard  deviation-to-mean  ratio  of  the  field  of  a 
vertical  dipole  and  4-element  dipole  arrays  in  the  forest  with  eejf  =  1.03  +  /0.036, 
canopy  height  H  =  20m,  tree  density  0.05/m2,  and  tree  trunk  height  of  h  =  15m. 


Ground(8g) 


Figure  1:  Geometry  of  a  forest  model  for  characterization  of  wave  propagation  in  a 
forest  environment. 


27 


Canopy-air  interface 


Figure  2:  Wave  propagation  mechanisms  contributing  to  the  mean  field  without  tree 
trunks  in  a  forest  environment. 


a  forest  environment. 
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X 


Figure  5:  Global  and  local  coordinate  systems  used  in  the  MoM  formulation  of  sparse 
scatterers. 


Figure  6:  Application  of  reciprocity  theorem  for  the  computation  of  scattered  field  of 
cylinder  from  a  nearby  dipole  embedded  in  a  dielectric  slab  above  ground  plane. 


30 


0  50  100  150  200  250  300  350  400  450  500 

Distance  (m) 

#  (b) 

Figure  7:  Relative  percentage  error  in  magnitude  (a)  and  phase  (b)  of  scattered 
electric  field  from  50  infinite  dielectric  cylinders  that  are  located  along  y-axis  with 
£  2m  separation  when  a  TE  plane  wave  is  obliquely  incident  with  6  =  60°  at  50  MHz. 
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Figure  8:  Mean  and  standard  deviation  of  scattered  electric  field  generated  by  scat¬ 
tered  in  each  20m  bin  with  TM  source(x  =  20,  y  =  25)  and  an  observation  point(x 
=  780,  y  =  25).  Three  different  frequencies  are  used,  10,  30,  and  50MHz  and  total  # 

number  of  scatterers  is  2000.  (a)  Magnitude  of  mean  (b)  Magnitude  of  standard 
deviation.  32 


-38.95 


Figure  9:  Comparison  of  the  mean  field  (a)  and  standard  deviation  (b)  of  Ez  of  the 
reduced  problem  that  keeps  600  scatterers  near  the  source  and  250  scatterers  near  the 
observation  point  with  those  of  the  complete  problem  that  contains  2000  scatterers, 
with  TM  line  source  excitation. 
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(C) 

Figure  10:  Comparison  of  the  standard  deviation  for  TE  ((a)  and  (b))  and  TM  ((c)) 
excitations  of  the  complete  and  reduced  problems  that  keeps  200  (TE)  or  150  (TM) 
scatterers  near  the  source  and  150  scatterers  near  the  observation  point  with  those  of 
the  complete  problem  that  contains  500  scatterers,  with  TE  line  source  excitation. 


Figure  11:  The  ratio  of  the  fields  of  a  dipole  with  and  without  a  tree  trunk  as  a 
function  of  dipole  height  (z)  and  azimuthal  angle  around  the  cylinder.  The  dipole  is 
vertical  and  10cm  away  from  the  cylinder  surface  and  observation  point  is  1km  away. 
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9 


Spacing  in  A/2  increments 


(a) 


Figure  12:  Spatial  correlation  of  magnitude  (a)  and  phase  (b)  of  field  of  a  vertical 
dipole  (Ez)  along  x-axis  and  y-axis. 
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FrequencyfMHz] 


Figure  13:  Magnitude  of  mean  field  as  function  of  frequency  in  the  absence  of  tree 
trunks  and  the  ground  plane.  RMS  height  is  2m  and  ee//  =  1.03  +  j0.0036  and  the 
transmitter  and  receiver  are  17m  and  15m  below  the  canopy-air  interface,  respectively. 


xIO'5 


Figure  14:  Impulse  response  with  Gaussian  pulse  excitation  on  the  z  directed  dipole. 
40  non-uniformly  spaced  sampling  points  are  used.  Dot  line  is  free-space  response 
multiplied  by  the  path- loss  evaluated  at  50MHz  (see  Table  2). 
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4i 


Appendix  c:  Software  Description 

The  document  provides  a  short  summary  of  different  programs  that  are  needed  for  simulation  of 
wave  propagation  in  a  forest  environment.  There  are  32different  programs  and  input  files.  That 
different  simulations  can  easily  be  carried  out  by  combining  appropriate  modules.  Of  32  files  17 
are  header  files  that  contain  different  related  mathematical  routines  and  are  called  automatically 
from  the  main  programs.  Table  1  includes  the  names  and  functions  of  these  header  files.  All 
header  files  are  specified  by  a  suffix  “h”. 

1  -  Description  of  Program 


Table  1:  Header  Files 


File  Name 


Description 


CG_solver.h  Contains  routines  for  CG  type  matrix  solver 


build_and_solve_ 

system.h 


calculus. h 


complex.h 


constant_em.h 


em  tool.h 


Contains  routines  for  building  and  solving  matrix  equation  for  approxi¬ 
mate  MoM. 


Contains  routines  for  Gauss  quadrature  integration 


Contains  routines  for  complex  arithmetic 


Contains  physical  and  mathematical  constants 


Contains  routines  for  calculation  of  Fresnel  reflection  coefficient 


file_handler.h  Contains  a  routine  for  opening  files 


incident_field.h  Contains  a  routine  for  the  calculation  of  lateral  wave  without  tree  trunks 


mesh  tooll.h 


mom_tool.h 


my_malloc.h 


near_field.h 


random_tool.h 


read_conf.h 


reflected  wave.h 


sparse_matrix.h 


specialfun.h 


Contains  a  routine  for  generating  tree  locations  and  discretizing  each 
cylinder 


Contains  a  routine  for  the  calculation  of  impedance  matrix. 


Contains  a  routine  using  “malloc”  function 


Contains  a  routine  for  the  calculation  of  near  field  from  finite  cylinder 


Contains  random  number  generator  routine 


Contains  a  routine  for  reading  input  files 


Contains  a  routine  for  the  direct  and  reflected  scattered  field  from  the 
ground  plane  using  near_field.h 


Contains  a  function  for  indexing  sparse  matrices 


Contains  special  functions  routines 


There  are  also  5  programs  that  handle  pre-  and  post-processing.  These  program  files  are  denoted 
by  a  suffix  “.c”  and  listed  in  Table  2. 


Table  2:  Pre-  and  Post-processing  Files 


File  Name 

Description 

mesh_s_corr.c 

Generates  50  different  realization  of  tree  trunk  positions  for  spatial 
correlation  simulation 

mesh_time.c 

Generates  a  single  realization  of  tree  trunk  position  for  time  domain 
simulation 

file_merger_s_corr.c 

Performs  the  calculation  of  output  results  (mean,  standard  deviation, 
and  correlation  coefficient)  for  spatial  correlation  simulation 

file_merger_time.c 

Merges  output  files  of  time_domain.c  (see  Table  3)  into  a  single  file 
used  by  time_post.c 

time_post.c 

Calculates  time  domain  response  from  the  frequency  domain  results 

There  are  3  main  programs  that  simulate  wave  propagation  in  a  forest  environment.  These  pro¬ 
grams  are  specified  by  a  suffix  “.c”  and  listed  in  Table3. 


Table  3:  Main  Programs 


File  Name 

Description 

array.c 

Calculates  field  at  one  observation  point  with  one  dipole,  or  array  excita¬ 
tion 

s_corr.c 

Calculates  spatial  correlation  along  x-axis  or  y-axis  with  one  dipole,  or 
array  excitation 

time_domain.c 

Calculates  spectral  domain  solution  for  time  domain  response  with  one 
dipole,  or  array  excitation 

There  are  3  different  input  files  that  are  needed  to  run  the  main  program  and  are  summarized  in 
Table  4. 


Table  4:  Input  Files 


File  Name 

Description 

array.conf 

Used  for  array.c  and  s_corr.c,  and  time_domain.c 

Table  4:  Input  Files 


File  Name 

Description 

array_info.conf 

Contains  locations,  orientation,  and  phase  difference  of  the  array  ele¬ 
ments 

time_domain.conf 

Used  for  time_domain.c 

There  are  4  m-files  that  plot  the  results  of  simulations.  Table  5  includes  names  and  their  functions. 


Table  5:  m-Files 


File  Name 

Description 

array_m.m 

Plots  the  resulting  electric  field  of  “array.c”  as  function  of  realization 

s_corr_m.m 

Plots  the  results  of  “s_corr.c”  using  “sampling_m.m” 

sampling_m.m 

Matlab  function  for  interpolation  using  sampling  theorem 

time_domain_m.m 

Plots  the  results  of  “time_domain.c”. 

2  -  File  Format  of  Input  Files 

1)  array.conf 

•  The  input  values  should  be  inserted  under  the  variable  name  in  the  order  shown  below. 


Table  6:  array.conf 


Variable  name 

Description 

Float  or 
Interger 

frequency: 

frequency  [Hz] 

Float 

obs_x,  obs_y, 
obs_z 

x,  y,  and  z  axis  of  Observation  point  [m] 

Float 

radius 

radius  of  cylinder  [m] 

Float 

min_gap 

minimum  distance  (min_gap*  wavelength)  between  cylinders 

Float 

density 

2 

density  of  cylinders  [lm  ] 

Float 

forest_H 

height  of  forest  [m] 

Float 

Table  6:  array.conf 


Variable  name 

Description 

Float  or 
Interger 

cyl_h 

Mean  height  of  cylinder  [m] 

Float 

cyl_v 

Variation  of  height  of  cylinder  [m] 

Float 

error_c 

error  criterion  for  iteration  solver 

Float 

er_eff_r 

Real  part  of  effective  dielectric  constant  of  the  canopy 

Float 

er_eff_i 

Image  part  of  effective  dielectric  constant  of  the  canopy 

Float 

eg_r 

Real  part  of  dielectric  constant  for  ground 

Float 

eg_i 

Image  part  of  dielectric  constant  for  ground 

Float 

eg_sigma 

Conductivity  of  dielectric  constant  for  ground 

Float 

#of_antenna 

Number  of  antennae  in  the  array  under  consideration 

Integer 

#of_mesh 

Number  of  discretization  points  along  a  cylinder  diameter 
used  in  the  method  of  moments. 

Integer 

#of_scattererl 

Number  of  scatterers  near  source 

Integer 

#of_scatterer2 

Number  of  scatterers  near  receiver 

Integer 

#of_iteration 

Number  of  iteration  number 

Integer 

#of_points 

Number  of  segments  along  z-axis  of  cylinders  for  numerical 
integration.  (>  L/0.15A,) 

Integer 

Maxit 

number  of  iteration  for  iterative  solver,  (typically  around  3 
times  the  matrix  size) 

Integer 

Note:  Float  type  variables  should  be  input  like  1.0  or  1.. 
2)array_info.conf 


Table  7:  array_info.conf 


Variable  name 

Description 

Float  or 
Interger 

ant_x,  ant_y,  ant_z 

Coordinate  of  the  transmitter  antenna  [m] 

Float 

ori_x,  ori_x,  ori_x 

Orientation  of  transmitter 

Float 

Table  7:  array_info.conf 


Variable  name 

Description 

Float  or 
Interger 

phase_diff 

Phase  difference  with  respect  to  the  first  antenna 
[Degrees]. 

Float 

Note:  The  coordinate,  orientation  and  phase-difference  for  each  antenna  in  the  array  should 
appear  in  separate  successive  lines. 

Example  -  The  input  file  for  three  z  directed  antennas  in  an  array  with  a  progressive  phase  180 
degrees. 


ant_x 

ant_y 

ant_z 

ori_x 

ori_y 

ori_z 

phase_dif 

20. 

25. 

3. 

0. 

0. 

1. 

0. 

20. 

30. 

3. 

0. 

0. 

1. 

180. 

20. 

35. 

3. 

0. 

0. 

1. 

360. 

3)  time_domain.conf 


Table  8:  time_domain.conf 


Variable  name 

Description 

Float  or 
Interger 

starting_f 

Starting  frequency  [Hz] 

Float 

end_f 

Stop  frequency  [Hz] 

Float 

sN 

Total  number  of  samples  between  starting_f  and  end_f 

Integer 

dis_.ro 

Distance  between  transmitter  and  receiver  [m] 

Float 

3  -  Running  the  Programs 

For  “array.c”,  simply  compile  and  run. 

For  “s_corr.c”  and  “time_domain.c”,  follow  the  steps  shown  below. 

1st  step:  To  generate  tree  locations  use  “mesh_s_corr.c”  for  spatial  correlation  calculation  and 
“mesh_time.c”  for  time  domain  calculation.  The  “mesh_s_corr.c”  generates  50  samples  of  tree 
arrangements  for  Monte  Carlo  simulation,  “meshjime.c”  generates  only  one  sample  of  tree 
arrangement.  The  “meshjime.c”  receives  input  argument  as  “meshjime  n”  where  n  means 
nth  realization. 


2nd  step:  Run  main  program,  “s_corr.c”  or  “time_domain.c”. 

“s_corr.c”  computes  fields  at  9  different  equally-spaced  points  along  a  line(specified  in 
“mesh_s_corr.c”).  Separation  between  two  adjacent  points  is  A/2  where  A,  is  wavelength. 

To  make  the  program  conducive  for  parallel  computation,  s_corr  is  run  using  three  parameters. 
To  run,  type  “s_corr  n  i  j”  with  n,  i,  j  being  integer  and  positive  numbers,  n  specifies  the  real¬ 
ization  number  ns  {0,  iV  -  1 } ,  i  specifies  the  observation  point  is  {1,9},  and  j  s  {1,2}  is 
a  parameter  specifying  computation  for  scatterers  near  the  source(j  =  0)  and  near  the  observa¬ 
tion  point(j  =  1).  To  compile,  type  “gcc  -o  s_corr  s_corr.c  -lm”,  which  generates  an  executable 
file  called  s_corr.  For  Monte-Carlo  simulation  the  following  procedure  may  be  followed: 

for  n  =  0:N-1 
for  i  =  1:9, 
for  j  =  1:2, 
s_corr  n i j 


where  N  is  the  total  number  of  realization. 

Note:  You  can  change  the  direction  of  observation  point  by  changing  “const  char  direction  = 
’x’;”  in  the  “mesh_s_corr.c”  that  means  observation  point  moves  along  x-axis.  If  you  change 
this  sentence  to  “const  char  direction  =  ’y’;”  this  means  observation  point  moves  along  y-axis. 

The  program,  “time_domain.c”  calculates  fields  at  sN  different  frequency  points  which  is  spec¬ 
ified  in  “time_domian.conf  ’  that  are  determined  by  the  Gaussian  quadrature  procedure.  To  run 
this  program,  type  “time_domain  n”  with  n  being  interger  and  positive  number. 
n  e  {0,  sN  -  1 }  is  a  index  pointing  a  frequency  point. 

3rd  step:  After  finishing  the  2nd  step,  run  “file_merger_s_corr.c”  for  spatial  correlation  or 
“file_merger_time.c  and  then  time_post.c”  for  time  domain  analysis  to  obtain  the  final  results. 

For  Monte-Carlo  simulation  the  procedure  shown  below  may  be  followed: 

for  n=  1:N, 
mesh_time  n 
for  i  =  0:39, 
time_domain  i 
end 

file_merger_time 

time_post 

move  output  files  to  other  file  name  to  prevent  overwriting  those  at  next  simulation 
end 


4  -  Output  Files 

1)  array.c 

There  are  3  output  files,  dipole_x_p_my.dat,  dipole_y_p_my.dat,  and  dipole_z_p_my.dat. 


dipole_x_p_my.dat  includes  the  x-component  of  scattered  electric  field 
dipole_y_p_my.dat  includes  the  y-component  of  scattered  electric  field 
dipole_z_p_my.dat  includes  the  z-component  of  scattered  electric  field 

Each  line  of  the  output  files  include  the 

real  part  and  the  imaginary  part  of  the  resulting  field  respectively. 


2)  s_corr.c 

There  are  3  output  files  for  this  program  as  well.  These  are  s_corr_mean.dat,  s_corr_var.dat, 
and  s_corr.dat. 

s_corr_mean:  Mean  field  at  the  observation  points 
s_corr_var.dat:  Standard  deviation  of  field  at  the  observation  points 
s_corr.dat:  Spatial  correlation  coefficient  with  respect  to  the  center  point 

The  format  of  data  in  files,  s_corr_mean.dat  and  s_corr.dat  is  as  follows. 

Re[Ex]  Im[Ex]  Re[Ey]  Im[Ey]  Re[Ez]  Im[Ez] 


The  format  of  data  in  file,  s_corr_var.dat  is  of  the  following  form. 
Var[Ex]  Var[Ey]  Var[Ez] 


Note:  The  program  “s_corr_m.m”  is  an  m-file  which  plots  spatial  correlation  after  interpolat¬ 
ing  data  using  the  sampling  theorem. 

3)  time_domain.c 

There  are  6  output  files  for  this  program. 
ifx_TD_inh.dat:  x-component  of  lateral  wave 
ify_TD_inh.dat:  y-component  of  lateral  wave 
ifz_TD_inh.dat:  z-component  of  lateral  wave 
sfx_TD_inh.dat:  x-component  of  scattered  wave 
sfy_TD_inh.dat:  y-component  of  scattered  wave 
sfz_TD_inh.dat:  z-component  of  scattered  wave 

The  format  of  data  in  these  files  are  as  follows: 

real  part  of  the  resulting  field  imaginary  part  of  the  resulting  field 


5)  Description  of  main  variables 


Table  9:  Main  Variables 


Variable 

Description 

Type 

N 

Total  number  of  elements  in  a  cylinder 

Integer 

Table  9:  Main  Variables 


Variable 

Description 

Type 

mN 

Total  size  of  matrix 

Integer 

N1 

Number  of  discretization  points  along  a  cylinder  diameter 
used  in  the  method  of  moments 

Integer 

i_N 

Number  of  iteration 

Integer 

S_N 

Total  number  of  scatterers  =  S_N1  +  S_N2 

Integer 

S_N1 

Number  of  scatterers  near  transmitter 

Integer 

S_N2 

Number  of  scatterers  near  receiver 

Integer 

maxiter 

Max.  number  of  iteration  for  iterative  solver 

Integer 

int_N 

Number  of  segments  along  z-axis  of  cylinders  for  numerical 
integration 

Integer 

a_N 

Number  of  antennae 

Integer 

kr 

Float 

kO 

Propagation  constant  in  free  space 

Float 

y_o 

Free  space  admittance 

Float 

zO 

Free  space  impedance 

Float 

lambda 

Wavelength 

Float 

ww 

Weight  factor  for  Gauss  quadrature 

Float  array 

XX 

Abscissas  point  for  Gauss  quadrature 

Float  array 

er_r 

Real  part  of  effective  dielectric  constant  of  the  canopy 

Float 

er_i 

Image  part  of  effective  dielectric  constant  of  the  canopy 

Float 

er_r 

Real  part  of  dielectric  constant  for  ground 

Float 

eg_i 

Image  part  of  dielectric  constant  for  ground 

Float 

sigma2 

Conductivity  of  dielectric  constant  for  ground 

Float 

fO 

Frequency 

Float 

wO 

2nf0 

Float 

j_error 

Error  criterion  for  iterative  solver 

Float 

ra 

Radius  of  cylinder 

Float 

Table  9:  Main  Variables 


Variable 


Description 


Type 


min_gap 


density 


x_c,  y_c, 
z_x 


cx,  cy,  cz 


array_posi 


array_ori 


array_pha 


height_s 


center 


delta_x, 

delta_y 


Z1 


im_field[3] 


pre_constl, 

pre_constxy 


minimum  distance  (min_gap*wavelength)  between  cylinders 

Float 

Density  of  cylinders 

Float 

Height  of  forest 

Float 

Mean  height  of  cylinder 

Float 

Variance  of  height  of  cylinders 

Float 

Coordinate  of  the  receiver 

Float 

Center  coordinate  of  each  cylinder 

Float  array 

Location  of  each  array  antenna 

Float  array 

Orientation  of  each  array  antenna 

Float  array 

Phase  difference  of  each  array  antenna 

Float  array 

Height  of  each  cylinder 

Float  array 

Coordinate  of  each  element 

Float  3-D  array 

Size  of  each  element  in  x-,  and  y-  axis 

Float 

Impedance  of  effective  dielectric 

Complex 

Scattered  field  from  scatterers  near  the  source. 

0:  x-comp. 

1:  y-comp. 

2:  z-comp. 

Complex  array 

Dielectric  constant  of  the  effective  media 

Complex 

Dielectric  constant  of  the  ground 

Complex 

Wave  number  of  the  effective  media 

Complex 

Wave  number  of  the  ground 

Complex 

kz 

Complex 

Pre-calculated  pre-factor  for  calculating  approximate  MoM 
matrix 

Complex 

pre_constxy 

1 


Table  9:  Main  Variables 


Variable 

Description 

Type 

mat 

Diagonal  block  matrix 

Complex  array 

F 

Current  Vector 

Complex  array 

J_dis 

Current  distribution,  resulting  vector 

Complex  array 

P,PP 

Needed  memory  for  iterative  solver 

Complex  array 
or  2-D  array 
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//  One  dipole  or  array  simulation  in  forest  enviroment  #endif 
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Complex ( temp, 0.);  fcomplex  *data;  double  *t;  double  t0; 
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